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We study the onset of patterns in vertically oscillated layers of frictionless dissipative particles. 
Using both numerical solutions of continuum equations to Navier-Stokes order and molecular dynam- 
ics (MD) simulations, we find that standing waves form stripe patterns above a critical acceleration 
of the cell. Changing the frequency of oscillation of the cell changes the wavelength of the resulting 
pattern; MD and continuum simulations both yield wavelengths in accord with previous experimen- 
tal results. The value of the critical acceleration for ordered standing waves is approximately 10% 
higher in molecular dynamics simulations than in the continuum simulations, and the amplitude 
of the waves differs significantly between the models. The delay in the onset of order in molecular 
dynamics simulations and the amplitude of noise below this onset are consistent with the presence 
of fluctuations which are absent in the continuum theory. The strength of the noise obtained by fit 
to Swift-Hohenberg theory is orders of magnitude larger than the thermal noise in fluid convection 
experiments, and is comparable to the noise found in experiments with oscillated granular layers 
and in recent fluid experiments on fluids near the critical point. Good agreement is found between 
the mean field value of onset from the Swift-Hohenberg fit and the onset in continuum simulations. 
Patterns are compared in cells oscillated at two different frequencies in MD; the layer with larger 
wavelength patterns has less noise than the layer with smaller wavelength patterns. 

PACS numbers: 45.70.Qj,05.40.Ca,47.54.+r 



I. INTRODUCTION 
A. Background 

A successful hydrodynamic theory of granular media 
could allow scientists and engineers to exploit the pow- 
erful techniques of fluid dynamics to describe granular 
phenomena. Recent experiments [l], [1J and simulations 
[3| demonstrate the potential for hydrodynamic theory 
to describe granular media; however, the validity of such 
methods has not yet been established for a general de- 
scription of granular flow phenomena [J, @ . 

Several proposed rapid granular flow models use equa- 
tions of motion for continuum fields - number density 
n, velocity u, and granular temperature T (|T is the 
average kinetic energy due to random particle motion) 
0> H Q • In one approach, particle interactions are mod- 
eled with binary, inelastic hard-sphere collision operators 
in kinetic theory to derive continuum equations to Euler 
[icj |. Navier-Stokes [HI], and Burnett [l2j order. In this 
paper, we use 3D simulations of continuum equations to 
Navier-Stokes order and 3D inelastic hard-sphere molec- 
ular dynamics (MD) simulations to investigate the onset 
of standing wave patterns in vertically oscillated granular 
layers. 



B. Standing wave patterns in oscillated granular 
layers 

Vertically oscillated layers have provided an important 
testbed for granular research. Flat layers of grains on 



a plate oscillating sinusoidally in the direction of grav- 
ity exhibit convection [l3| , clustering (bj, shocks [l5[ , 
steady-state flow fields far from the plate [n| , and stand- 
ing wave pattern formation [TtJ ■ 

A layer of grains on a plate oscillating sinusoidally in 
the direction of gravity with frequency / and amplitude 
A will leave the plate at some time in the cycle if the 
maximum acceleration of the plate is greater than that 
of gravity. The layer dilates above the plate, then col- 
lides with the plate later in the cycle and is compressed 
on the plate by this collision. Above a critical value of 
acceleration, standing wave patterns spontaneously form 
in the layer. This pattern is subharmonic with respect to 
the plate, repeating every 2/f [P?! ]. 

Various subharmonic standing wave patterns, includ- 
ing stripe, square, and hexagonal patterns, have been 
found experimenta lly, d epending on the nondimensional 
frequency /* = / \J H j g and the nondimensional acceler- 
ational amplitude V = A (27r/) 2 /<?, where H is the depth 
of the layer as poured, and g is the acceleration due to 
gravity [17] . 

Studies using hydrodynamic equations have not yet 
yielded the standing wave patterns observed in experi- 
ments. Here we investigate the onset of ordered standing 
wave patterns using fully three-dimensional (3D) simula- 
tions of continuum equations to Navier-Stokes order as 
well as molecular dynamics (MD) simulations. We use a 
continuum model for frictionless, inelastic particles, and 
investigate the onset of stripe patterns. 
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C. Fluctuating hydrodynamics 

Near the onset of convection patterns in Raylcigh- 
Benard convection of fluids, fluctuations caused by ther- 
mal noise create deviations from dynamics predicted from 
linear theory. These fluctuations are described by the ad- 
dition of terms to the Navier-Stokes equations; this the- 
ory is known as fluctuating hydrodynamics [3, [lfl H(| . 
Recent experiments have shown that fluctuating hydro- 
dynamics theory accurately describes the dynamics of flu- 
ids near the onset of convection [22|, [23[ • 

Experimental investigations of coherent fluctuations 
and pattern formation in oscillated granular layers have 
indicated that fluctuations due to the movement of in- 
dividual grains play a much more significant role in the 
collective behavior of granular media than do thermal 
fluctuations in ordinary fluids 24j. Thus, a consistent 
theory of granular hydrodynamics may need to include 
fluctuations. 



D. Model system 

We simulate a layer of grains on an impenetrable plate 
which oscillates sinusoidally in the direction of gravity. 
The layer depth at rest is approximately H — 5.4er, where 
the grains are modeled as identical, frictionless spheres 
with diameter a and coefficient of restitution e = 0.7. 
For most of the paper, we study the onset of patterns as 
a function of T, while the frequency of plate oscillation 
is held constant at /* = 0.4174. This corresponds to a 
frequency of 56 Hz for particles with a diameter of 0.1 
mm. For T > 2.5, stripes are seen experimentally for 
a range of parameters, including /* = 0.4174, H = 5.4 
fl7| . In Sec. IIIIBI and Sec. IIVCI frequency is varied to 
investigate the effect of changing frequency on pattern 
formation. 

Experiments (2f| and MD simulations [26] indicate 
that inter-particle friction plays an important role in the 
standing wave patterns. MD simulations with friction be- 
tween particles have quantitatively reproduced the stripe, 
square, and hexagonal subharmonic standing wave pat- 
terns seen experimentally for a wide range of parameters 
[27| . However, MD simulations using frictionless parti- 
cles do not yield stable square or hexagonal patterns, but 
only yield stripe patterns, and exhibit the onset of pat- 
terns at lower F than that seen for frictional particles 
[26| . This result is consistent with experiments which 
show that reducing friction by adding graphite can de- 
stabilize square patterns [25] . In this study, we neglect 
the effects of friction in our continuum and MD simu- 
lations, and study only the onset of stripe patterns in 
frictionless layers. To investigate other patterns such as 
squares or hexagons, simulations would need to include 
friction between particles. 

We use MD and continuum simulations to investigate 
the dynamics of this system near onset, and use simu- 
lations of the Swift-Hohenberg (SH) model equation to 



compare our results between the two. Section II describes 
the methods used to simulate and analyze these patterns, 
Sec. Ill compares patterns formed in MD and contin- 
uum simulations. Section IV compares MD simulations 
to Swift-Hohenberg theory, and Sec. V presents our con- 
clusions. 



II. METHODS 
A. Molecular dynamics simulation 

We use an inelastic hard sphere molecular dynamics 
simulation, which was previously used in conjunction 
with the continuum simulation used in this paper to 
model shock waves in a granular shaker (28j . This same 
MD code with friction added has been found to describe 
well the patterns observed in experiments on oscillating 
granular layers [27l[29j. 

The collision model assumes instantaneous binary col- 
lisions in which energy is dissipated, as characterized by 
the normal coefficient of restitution e. We neglect sur- 
face friction between particles, as well as between the 
particles and the plate. To prevent inelastic collapse, 
we use a coefficient of restitution which depends on the 
relative colliding velocity of the particles v n : e(v n ) — 

1 — 0.3 (y n / \/ga) 3 ^ 4 for v n < yfgcf, and e = 0.7 otherwise 

The MD simulations are calculated in a box of size 
L x x L y in the horizontal directions x and y, where L x and 
L y are varied to investigate patterns with different wave- 
lengths. The simulations use periodic boundary condi- 
tions in the horizontal directions, an impenetrable lower 
plate which oscillates sinusoidally between z = and 
z = 2A, and an upper plate fixed at a height z = 200c, 
as in the previous investigation of shock propagation [28| . 

For each MD simulation, (L x /<j) x (L y /a) x 6 particles 
were used. In actual packings seen experimentally, 6/cr 2 
particles per unit area of the bottom plate corresponds 
to a layer depth H — 5Aa as poured, representing a 
volume fraction v « 0.58. [27]. The total mass of the 
layer matches that of the continuum simulations. 



B. Continuum simulation 

We use a continuum simulation previously used to 
model shock waves in a granular shaker [28[. Our sim- 
ulation numerically integrates continuum equations of 
Navier-Stokes order proposed by Jenkins and Richman 
jll] for a dense gas composed of frictionless (smooth), 
inelastic hard spheres. We integrate these hydrodynamic 
equations to find number density, momentum, and gran- 
ular temperature, using a second order finite difference 
scheme on a uniform grid in 3D with first order adaptive 
time stepping [28[. 

As in our MD simulations, the granular fluid in the 
continuum simulations is contained between two impen- 
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etrable horizontal plates at the top and bottom of the 
container, where the lower plate oscillates sinusoidally 
between height z = and z = 2A. In our MD simula- 
tions, the ceiling is fixed in space at a height of z = 200<7, 
but to minimize computation time, the ceiling in con- 
tinuum simulations is located at height 80a above the 
lower plate and oscillates with the bottom plate. In our 
previous study of shock formation, changing the ceiling 
height from 200a to 80<r resulted in a fractional root 
mean square difference of less than 1% in the shock lo- 
cation over one cycle [28j . 

As in the MD simulations, we use periodic horizontal 
boundary conditions and boxes of size L x x L y in the 
horizontal directions x and y, where L x and L y are varied. 
In each case, continuum simulations are compared to MD 
simulations with the same horizontal dimensions L x and 
L y . The numerical methods, boundary conditions at the 
top and bottom plate, and grid spacing are the same as 
used in the previous study of shocks [281 ] . 

The energy loss due to collisions in continuum simula- 
tions is characterized by a single parameter, the normal 
coefficient of restitution e = 0.70. Throughout this pa- 
per, we use units such that particles in MD simulations 
have mass unity, and the total mass of the layer in the 
continuum simulations matches that used in MD simula- 
tions. 



C. Characterizing patterns 

To visualize peaks and valleys formed by standing wave 
patterns, we calculate the height of the center of mass of 
the layer, z cm (x, y, t) as a function of horizontal location 
in the cell at various times in the cycle. At a given time 
to and horizontal location (xo,yo), z cm (xo,yo,to) is the 
center of mass of all particles whose horizontal coordi- 
nates lie within a bin of size Ax i n x Ay^m centered at 
(^Oj^o)- F° r continuum simulations, we use the simula- 
tion grid size to define the bins: Axbin = Ax = 2a and 
Aybm = Ay = 2a. For MD simulations, we use bins 
of size 2a x 2a in Section IIIII to compare to continuum 
simulations with the same bin size. Peaks in the pattern 
correspond to maxima of z cmi and valleys correspond to 
minima. 

To measure the amplitude of patterns and fluctuations 
in continuum and MD simulations, we examine the devi- 
ation of the height of the center of mass of the layer as a 
function of horizontal location in the cell from the center 
of mass height averaged over the cell at that phase in the 
cycle: 



tp(x, V, t) = z cm (x, y, t) - (z cm (x, y, t)) 



(1) 



where x and y are the horizontal coordinates, t is the time 
in the cycle, z cm (x, y) is the height of the center of mass 
of the layer at horizontal location (x, y), and the brackets 
represent an average over all horizontal locations in the 
cell at a given time t. 



Throughout this paper, we characterize the patterns 
at the beginning of a sinusoidal oscillation cycle, such 
that the plate is at its equilibrium position and moving 
upwards. Using this definition, (?/? 2 (i)) represents the 
mean square deviation of the height of the layer from the 
mean height of the layer at that phase of the plate. We 
note that (ip 2 ) is large for layers with high amplitude 
patterns or fluctuations, and goes to zero as the layer 
becomes perfectly flat. 

In addition to (ip 2 ), we distinguish between or- 
dered patterns (stripes) and disordered fluctuations by 
characterizing the long range order of the pattern. 
To characterize the long range order of the patterns, 
we first calculate the power spectrum of the pattern: 



S (k x , ky, t) 



tp (k x , ky, t) 



Vy , t) 



where ip (k x , k„ 

We then transform 

— / £-2 i L.2 
r — \ I ru*. l 'Vy ; 



Jq X J q " ip (x, y, t) e~ tkxX e~ ik v y dxdy 

to polar coordinates in k space: k r = ^J~k% + k^, kg — 

tan^ 1 (k y /k x ) to find S(k r , kg) in the range < kg < n. 
We integrate radially to find the angular orientation of 
the power spectrum: S(kg) = J S (k r ,kg) k r dk r , where 
K = 2?r ^ t "" . We bin kg into 21 bins between kg = and 
kg = 7T, and characterize the long range order of the pat- 
terns by the fraction of the total integrated power that 
lies in the bin with the maximum power: 



Pn 



(0) 



JZS^dkg' 



(2) 



where S max (9) is the integrated power within an angle 
7r/21 of the maximum value of S(9). Thus P m ax is the 
fraction of the total power that lies within approximately 
7r/21 of the angular location of the maximum power. For 
a perfectly disordered state, with equal power in all direc- 
tions, P m ax would approach ^- ss 0.05, while P max = 1 
for a state with all power in a single bin. Thus P max 
provides a measure of order when stripes form. 



III. PATTERN ONSET AND DISPERSION 
A. Stripe patterns 

Experimental investigations of shaken granular layers 
have shown that above a critical acceleration of the plate 
r c , standing wave patterns form spontaneously. These 
patterns oscillate subharmonically, repeating every 2/f, 
so that the location of a peak of the pattern becomes a 
valley after one cycle of the plate, and vice versa [l7T |. 

Continuum and MD simulations produce standing 
wave patterns for T = 2.2 and /* = 0.4174 (Fig.rrj. Al- 
ternating peaks and valleys form a stripe pattern which 
oscillates at f /2 with respect to the plate oscillation; a 
location in the cell which represents a peak during one 
cycle will become a valley the next cycle, and then re- 
turn to a peak on the following cycle. For a box of size 
126(7 x 126cr in the horizontal direction, six wavelengths 
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FIG. 1: An overhead view of a layer of grains, showing the 
center of mass height z cm as a function of horizontal posi- 
tion (x,y) in a cell with horizontal dimensions L x x L y = 
126<t x 126<t, from (a) MD simulations and (b) continuum 
simulations. Peaks of the layer corresponding to large center 
of mass height z cm are shown in white; valleys corresponding 
to low Zcm are shown in black. 
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FIG. 2: Dispersion relation for stripes which form perpendic- 
ular to the long dimension of cells with horizontal dimensions 
168<r x lOcr. Data for continuum simulations are shown as 
triangles and MD simulations as circles; points where con- 
tinuum and MD simulations yield the same wavelength are 
shown as squares. In both continuum and MD simulations, 
the dominant wavelength of the final oscillatory state A fits 
very well to the dispersion relation found in experiments 
A* = 1.0 + i.i/-i-3 2 ±0 03 ( solid line ) r£ 2 j. Error bars in both 
simulations are calculated exclusively from discretization due 
to periodic boundary conditions in a finite size box. 



lit in the box in both MD and continuum simulations, 
yielding a wavelength of 21a ± 4cr in both continuum and 
MD simulations (Fig. [TJ . 



B. Dispersion Relations in Continuum, MD, and 
Experiment 

Experiments have shown that the wavelength A of 
standing wave patterns in shaken granular layers depends 
on the frequency of the plate oscillation [3CJ, [3l|, |32| . 
For a range of layer depths and oscillation frequencies, 
experimental data for frictional particles near the on- 
set of patterns were found to be fit by the function 
A* = 1.0 + i.i/*-i-32±o.03 j where A * = X J H 

We investigate the frequency dependence of standing 
waves in continuum simulations and in MD simulations 
of frictionless particles. Dimensionless accelerational am- 
plitude r = 2.2 was held constant while dimensionless 
frequency /* was varied. Simulations were conducted in 
a box of horizontal extent L x — 168(7 and L y = lOcr. This 
orientation causes stripe patterns to form parallel to the 
y— axis. The dominant wavelength in each case was cal- 
culated from S (k x , k y ,t) by finding the wavenumber k x 
in the x— direction which exhibited the maximum power 
during one cycle of the oscillatory state of the pattern. 
Due to the periodic boundary conditions and finite box 
size, wavelengths must fit in the box an integer number 



of times. This finite size effect of quantized wavelength 
yields inherent uncertainty in the wavelength that would 
be selected in an infinite box. 

Wavelengths found in continuum and MD simulations 
are compared to the dispersion relation fit to experimen- 
tal data in Fig. [2] Investigation is limited to /* > 0.15 
by the box size, as only two wavelengths fit in the box in 
continuum simulations at this frequency. Neither simu- 
lation produced patterns for this box size for /* > 0.45. 
Both simulations agree quite well with the experimental 
fit throughout the range 0.15 < /* < 0.45. 

Comparison to the experimental fit shows that both 
MD and continnum simulations produce wavelengths 
consistent with experimental results for frictional parti- 
cles. These data indicate that friction seems to be unim- 
portant in wavelength selection through this parameter 
range. 



C. Layers Above and Below the Onset of Patterns 

Continuum and MD simulations exhibit pattern forma- 
tion above a critical acceleration of the plate; however, 
standing wave patterns are not observed below a critical 
value of T (Fig. |3J>. For T = 2.2, both MD (Fig. [3(b)] ) 



and continuum (Fig 
peaks and valleys w 



3(d)[ ) simulations show well defined 



lich form stripe patterns with two 



wavelengths fitting in the box of size L x 



42cr. 
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FIG. 3: An overhead view of the layer of grains, showing the 
center of mass height z cm (x,y) of the layer as a function of 
location in the box, for (a) MD simulations with a plate accel- 
eration with respect to gravity F = 1.9, (b) MD simulations 
with F = 2.2, (c) continuum simulations with F = 1.9, and 
(d) continuum simulations with V = 2.2. Peaks corresponding 
to large z C m are shown in white, while valleys corresponding 
to small Zcm are shown in black. The grayscale for all four 
images is given on the right. 




FIG. 4: The mean square deviation (^ 2 ) of the local center 
of mass height from the average center of mass height of the 
entire layer as a function of accelerational amplitude F for 
MD (triangles) and continuum (circles) simulations. The ver- 
tical dotted line represents the onset of stripe patterns in the 
continuum simulations. 



D. Onset of patterns in continuum simulations 



The only difference between this system and that investi- 
gated in Sec. IIII Al is the horizontal size of the cell; these 
patterns look very similar to a section of the patterns 
formed in the larger cell (Fig. [1]). Reducing the acceler- 
ational amplitude to T = 1.9 while keeping all other pa- 
rameters constant yields no ordered waves in either MD 



(Fig. 3(a)) or continuum (Fig. 3(c)). Thus both contin- 



uum and MD simulations appear to have a critical value 
of r somewhere in the range 1.9 <T C < 2.2, such that no 
patterns are formed for T < T c , and patterns are formed 
for r > r c . This critical value is lower than that found 
in experiments with frictional particles, where a similar 
onset of patterns is found at a critical value of T w 2.5 

Despite the similarities, differences between MD and 
continuum simulations are observable. For T = 1.9, 
the continuum simulation yields a very smooth, flat 
layer (Fig. 3(c)[ ), while MD exhibits visible fluctuations 
(Fig. 3(a)). For T = 2.2, the continuum simulations pro- 
duce stripes (Fig. |3(d)[ ) which are much smoother than 
those found in MD simulation (Fig. |3(b)[ ) . 

To explore the differences between the two simulations, 
we investigate the onset of patterns in more detail in 
continuum simulations and MD simulations separately. 



We investigate the onset of patterns in continuum sim- 
ulations by determining (tp 2 ) of standing waves for differ- 
ent values of T. Each simulation begins with a flat layer 
above the plate with small amplitude random fluctua- 
tions. The simulation is run until it reaches a periodic 
state, at which point (if) 2 ) is calculated as an average 
over ten cycles of the same phase of the plate. 

For r < 1.95, the initial fluctuations decay rapidly 
until the layer is quite flat, as represented by negligible 
values of (if> 2 ) (Fig- HI). As T increases, there is a sudden 
onset to large amplitude waves, as seen by the sudden 
jump in (ip 2 ) in Fig. |4j This onset occurs at the critical 
value T c = 1.955 ±0.005. For T < T c , initial fluctuations 
decay until the layer is very flat, while for all layers above 
onset (r > r c ), these waves produce ordered patterns of 
stripes similar to those in Fig. |3(d)| 



E. Onset of patterns in molecular dynamics 
simulations 



We examine the onset of patterns in MD simulations 
using the same methods as for the continuum equations. 
Figure 0] shows the mean square height deviation (fp 2 ) as 
a function of T for MD simulations as well as for contin- 
uum simulations. For each value of T, the simulation was 
run for 400 cycles of the plate until the layer reached a 
periodic state, then (^ 2 ) and P ma x were calculated from 
an average of the next 100 cycles. 

As in continuum simulations, ("0 2 ) grows with in- 
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creasing V. Unlike the continuum results, ("0 2 ) is non- 
negligible in MD simulations even for T < 1.95. There is 
still a sharp increase in the slope of the curve, but it is 
delayed until T > 2.1. 



IV. THE ROLE OF FLUCTUATIONS 

The MD simulations display an onset of ordered stripes 
that is delayed with respect to those found in continuum, 
and exhibit non-negligible (ip 2 ) even below the onset of 
ordered stripes. Since the hydrodynamic model used in 
the continuum simulations does not include a stochastic 
noise term characteristic of fluctuating hydrodynamics, 
the differences between the continuum and MD simula- 
tions may be consistent with the presence of noise in the 
MD simulations due to the small number of particles per 
wavelength. To test the hypothesis that these differences 
are consistent with the presence of fluctuations in molec- 
ular dynamics simulations, we compare MD simulations 
to results from the Swift-Hohenberg model. 



A. Swift Hohenberg simulation 

The Swift-Hohenberg (SH) model was developed to de- 
scribe thermal noise-driven phenomena near the onset of 
long range order in Rayleigh-Benard convection [20l | . Re- 
cent experimental evidence suggests similar phenomena 
in shaken granular experiments can be interpreted using 
the methods of fluctuating hydrodynamics |24| . 

The SH model describes the time evolution of a scalar 
field VsH(x,t): 



3^J 



SH 



dt 



= e- 1 



SH 



4%H+V(x,t), (3) 



where e is the bifurcation parameter, and 77 is a stochastic 
noise term of strength F, such that (rj (x, t) rj (x', t')) = 
2F6 (x — x') 6 [t — t'). In the absence of stochastic noise 
(F = 0), called the mean field (MF) approximation, there 
is a sharp onset of stripe patterns with long range order at 
e = e f F = [23,113. For F ^ 0, the effect of noise is to 
delay the onset of long range (LR) order to a new critical 
value: e FR > 0. The delay in onset is characterized by 
Ae c = e^ R — e^ IF . In addition, the presence of noise 
creates fluctuations below the onset of long range order 
(e<e^). 

The Swift-Hohenberg simulation displays a forward bi- 
furcation to stripes at onset, while MD simulations show 
slight (< 1%) hysteresis [24|. A more complicated SH 
model [13] yields square patterns with hysteresis; how- 
ever, in this work we compare stripe formation in MD 
simulations a simpler model of the effects of noise near a 
bifurcation (Eq. [3]). 

We numerically solve the SH equation using the scheme 
described in (35[, with the number of gridpoints N = 
42 x 42, and periodic boundary conditions. We use in- 
tegration timesteps of 0.5, and the size of each gridspace 



in the horizontal directions Ax = Ay = 0.29 so that two 
wavelengths of the resulting pattern fit in the box, to 
match MD and continuum simulations. The simulation 
was allowed to run for 8,000 timesteps to reach a final 
pattern; then (V'sff) an< ^ Pmax were calculated from an 
average of the next 2,000 timesteps, in the same way as 
(ip 2 ) and Pmax were calculated for MD and continuum 
simulations in Section Hi CI 



B. Comparing Swift-Hohenberg and molecular 
dynamics simulations 

To find the strength of the noise and the mean field 
onset, we fit the SH model to the data from MD simula- 
tions (Fig. [5]) by varying three parameters: F, Ae c , and 
an overall scale factor, as in [23l. [24l|. 

Of the three parameters, only the noise strength F 
changes the overall shape of the curve. For a given F, 
the SH simulation is run for a range of —0.2 < e < 0.2; 
ipSH and P m ax are calculated from the steady state so- 
lution for each value of e and compared to MD simula- 
tions. For consistency, (ip 2 ) and P ma x are calculated for 
MD simulations from bins of size Axbin = Aybi n = a 
throughout this section, so that the number of bins in 
both MD and SH simulations is 42 x 42. Increasing the 
bin size to Axun = Ay^ = 2a does not change any of 
the fit parameters to within our uncertainty. 

Note (V'iff) m SH simulations is found as a function 
of control parameter —0.2 < e$H < 0.2, while in MD 
simulations, (^md) is found as a function of control pa- 
rameter 1.7 < r < 2.3. To compare the onset of the 
SH model to the onset in MD simulations, we define 



_ (r ~nMF\ /-pMF 
(MD - (1 - l c J/l c , 



where Tq F is the mean field 



onset of patterns, comparable to esn = 0. However, we 
do not know a priori the value of Tq . 

We find that (tp 2 ) changes relatively smoothly in MD 
and SH simulations, making it difficult to pinpoint an 
onset of patterns from (V' 2 ) alone. However, there is a 
distinct onset of long range order in the system (Fig. [5} . 
For low r in MD, the fluctuations are disordered (cf 



Fig. 3(a) ), while for higher T, standing wave stripe pat- 
terns are observed (c/Fig. [3(b)[ ). A clear transition from 
disordered fluctuations to an ordered stripe pattern is 
demonstrated by the sharp increase in P max as T crosses 
the critical value for long range order, determined from 
Fig. |5(b)| as T^ R = 2.15 ± 0.01. A similar tr ansitio n to 
ordered stripes is seen in SH simulations (Fig. |5(b)[ ). 

The onset of long range order is used to establish a 
correspondence between T and e. For MD simulations, 
we measure the onset of long range order as the point of 
sharpest increase in P max (Fig. 5(b) ). In SH simulations, 
Ae c represents the onset of long range order. We match 
the single point of steepest increase of P m ax between the 
two curves. The measured value Ae c in SH then predicts 
the mean field onset Tq F corresponding to e = 0. 

Once the relationship between T and e is determined, 
the overall scale factor for a given F is found by a least 
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FIG. 5: Comparison of MD simulations (triangles) to the 
Swift-Hohenberg model (solid lines) for (a) (ip 2 ), and (b) 
global ordering P m ax (Eq. [2}, as a function of control pa- 
rameter e (bottom axis) for SH, and Y (top axis) for MD. 
The parameters for SH simulations are noise strength F — 
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FIG. 6: The squared residual R 2 between (ipMr>) an d (wsh) 
as a function of the noise strength F used in SH simulations. 
The best least squares fit is given by F — (1.2 ± 0.2) x 10~ 2 . 



simulations, and as a function of T for MD simulations. 
The fit shows good agreement in (tp 2 ) below e = 
(Fig. Although the parameters are fit only in the 
range 1.7 < T < , agreement is reasonable even for 

The three parameter fit not only allows for agreement 
in (i/> 2 ), but also matches the measure of order P m ax in 



the SH model to that found in MD simulation (Fig. 5(b) ). 



In both MD and SH simulation, below the critical value 
of long range order, the fluctuations are disordered, lead- 
ing to a small value in P max . When T crosses the criti- 
cal value, P max jumps up significantly, and the observed 
patterns are ordered stripes. Below the onset of stripes, 
when the fluctuations are constantly shifting and chang- 
ing, there is significant uncertainty in finding the value 



(1.2 ± 0.2) x 10 and a delayed onset of long range or- r r, u ax. • j-i i i a u 

L. ,lr ' n nnA rT-n ,,J , ■ „ ■ ° , JV . of Pmax, as seen by the noisy curve on the plot. Above 



der e c = 0.094. The global ordering jumps sharply at 
e^ R = 0.094, corresponding to T% R = 2.15 in MD (the 
vertical dotted line in the figure), representing a transition 
to stripe patterns, while (ip 2 ) increases smoothly through 
that transition. This fit predicts a mean field onset value of 
rf F = 1.965 ±0.007, corresponding to e^ IF = (the vertical 
dashed line in the figure). 



this onset, however, the standing waves produce stable 
stripes, and P m ax plateaus and remains quite constant, 
with good agreement between MD and SH simulations. 
Finally, the mean field onset Tf F = 1.965 ± 0.007 pre- 
dicted by this fit agrees remarkably well with the critical 
value r c = 1.955 ± 0.005 found in our simulations of 
Navier-Stokes order continuum equations. 



squares fit between (V'lff) an d (V'md) f° r the range 
1.7 < T < T^ R (see Fig. |5(b)[ ). This minimization proce- 
dure gives the best possible fit for a given value of F. 
This entire procedure is repeated for vary- 
F, minimizing the squared residual R 2 = 



tag 



S {{^ 2 md) ~ (V'li?)) /N, where N is the number 
of bins (Fig. [S]). The best fit yields an onset of long 
range order at Ae c =0.94, corresponding to T^ R — 2.15. 



Figure 5(a) shows < ip > as a function of e for SH 



C. Effect of changing wavelength on strength of 



If the noise effects arise from the finite particle number 
in MD, we may expect that this effect will decrease in 
systems in which there are more particles per wavelength 
of pattern. Since the number of particles in a volume A 3 
increases with increasing wavelength, we investigate the 
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FIG. 7: Comparison growth of (ipsu) normalized by the 
mean center of mass height of the layer a 2 (ip 2 ) / iz cm ) 2 = 
((zcm — {zcm)) 2 ) / (z C m) 2 for MD simulations with /* = 0.25 
(squares) and /* = 0.4174 (triangles). The lower frequency 
displays much smaller fluctuations below the onset of patterns 
than does the higher frequency. 



effect of changing frequency on the onset of patterns in 
MD simulations. For cells of horizontal extent 168ctx 10a, 
layers shaken with a frequency /* = 0.25 form peaks 
with a dominant wavelength A = 42er, which is twice 
the wavelength found for patterns investigated at /* = 
0.4174 (see Fig.©. 

We examine layers shaken at /* = 0.25 in cells of size 
L x = L y = 2A = 84cr, while holding constant layer depth 
H = 5.4 and restitution coefficient e = 0.70. We vary 
r through the same range 1.7 < T < 2.3 investigated 
for /* = 0.4174 earlier in this paper. Figure [7] shows 
the growth of (tpsH/ normalized by the mean center 
of mass height of the layer squared er 2 (ip- 2 ) / (z C m) 2 = 



{zcm — (z C m)) ) / (z C m) for MD simulations with fre- 
quencies /* = 0.25 and /* = 0.4174. 

The lower frequency (/* = 0.25) exhibits a much 
sharper jump in (ipg H ) than that seen at /* = 0.4174. 
Below this onset, the curve is much flatter for /* = 0.25, 
while at /* = 0.4174, the curve increases much more 
gradually through onset. Proportionally smaller fluctu- 
ations compared to pattern amplitude is consistent with 
lower noise strength for /* = 0.25 than that found for 
/* = 0.4174. In addition, the rapid growth of peaks and 
valleys occurs at a smaller value of T for /* = 0.25, cor- 
responding to an onset even below the mean field onset 
T^f F for the larger frequency. 

We follow the same procedure as for /* = 0.4174 to 
fit the data from MD simulation to the Swift-Hohenberg 
model. We note that for frictional particles, square pat- 
terns are formed for /* = 0.25; in the absence of fric- 



FIG. 8: An overhead view of the layer of grains from MD 
simulations at /* = 0.25, for V = 1.7 and F = 2.2. Note how 
much less noise there is below onset here (F = 1.7) compared 
to the results for /* = 0.4174 in Fig. [3] The images show 
the center of mass height z cm (x,y) of the layer as a func- 
tion of location in the box. These MD simulations use a cell 
which is L x = L y = 84a in the horizontal directions. Peaks 
corresponding to large z cm are shown in white, while valleys 
corresponding to small z cm are shown in black. The grayscale 
for both images is given on the right. 



tion, peaks and valleys remain disordered, and no reg- 
ular square lattice forms in experiments or MD simula- 
tions [2a, Ha] (see Fig. [8j) . Thus the onset of long range 
order is ill defined in this case. However, this lower fre- 
quency exhibits a much sharper onset in the growth of 
(V'sh)) which is used to find Ae c . The best fit yields a 
noise term F = (4 ± 3) x 10~ 4 , and a mean field onset 
of = 1.85 ± 0.01. Our hydrodynamic simulations 

find the flat layer becomes unstable at r c = 1.84 ± 0.01, 
which again agrees well with the mean field onset found 
from the fit to the SH model. 

The noise strength at /* = 0.4174 is approximately 
30 times larger than the noise strength at /* = 0.25. 
This leads to qualitatively different behavior of (iPsh) 
near onset, yielding a smoother curve for the higher fre- 
quency and a sharper onset for lower frequency. Finally, 
the onset is barely delayed for the lower frequency, with 
Ae c = 0.01 for /* = 0.25, as compared to Ae c = 0.10 for 
/* = 0.4174. 

Thus a change in frequency which increases the wave- 
length at onset by a factor of 2 decreases the amount of 
noise by a factor of 30. For Raylcigh-Benard convection 
in ordinary fluids, the functional dependence of F on n, 
u, T, and A is known [3(|[37|]. However, this dependence 
is not known for oscillated granular layers. Future inves- 
tigation of the dependence of F on shaking parameters 
/*, r, and H, or on hydrodynamic variables n, u, T in 
experiment and MD simulations may provide informa- 
tion on the dependence of the noise strength F that can 
be used in continuum simulations. 



V. CONCLUSIONS 

We have shown that continuum simulations can de- 
scribe important aspects of pattern formation in granular 
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materials. For a nondimensional frequency /* = 0.4174, 
both MD and continuum simulations of granular materi- 
als form stripe patterns of the same wavelength above a 
critical value T > T c , and display no stripes for r < T c . 
Further, the two simulations yield the same dependence 
of wavelength on frequency. These wavelengths agree 
with the dispersion relation found experimentally for fric- 
tional particles. 

The effect of fluctuations has been examined in sim- 
ulations of the Swift-Hohcnberg model. The results de- 
duced for the mean field onset in MD simulations agree 
well with the actual onset in continuum simulations for 
both /* = 0.4174 and /* = 0.25. 

We find the strength of the noise to be F = 
(1.2 ±0.2) x 10" 2 for stripes at /* = 0.4174, and F = 
(4 ± 3) x 10~ 4 for disordered squares at /* = 0.25. The 
value determined in an experiment for a slightly shal- 
lower granular layer at /* = 0.28 was F = 3.5 x 10~ 3 
[24| . which is within the range of noise values obtained 
in this investigation. The smallest noise strength found 
for our granular system is comparable to the largest noise 
strength found thus far in experiments in ordinary fluids, 
which obtained F — 7.1 x 10~ 4 in measurements near 
the critical point, while values typical for convection are 
closer to 10~ 9 [23|]. For /* = 0.4174, the noise is strong 
enough to delay onset of long range patterns by 10% in 
MD simulation, and influences strongly the behavior of 
the system even more than 20% below this onset. Thus 
noise plays an important role in granular media near the 
onset of patterns. 



This study indicates that hydrodynamic theory holds 
promise for investigating and understanding pattern for- 
mation in granular flows. However, quantitative com- 
parisons between continuum theory and experiment will 
require the addition of noise terms into the equations. 
The addition of noise would be an important step to- 
wards using the powerful tools of hydrodynamic theory 
to investigate problems of pattern formation in granular 
materials. 

The absence of friction in these simulations restricts 
our investigation to stripe patterns. Simulations without 
friction have not yielded the square and hexagonal pat- 
terns seen in experiments with frictional particles [26|. 
Further research into pattern formation using continuum 
simulations should investigate the most effective way to 
incorporate friction between particles into the continuum 
simulations and should examine how the strength of fric- 
tion in the simulation affects pattern formation in the 
system. 
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